Coarse-Grain Model for Lipid Bilayer Self- Assembly and Dynamics: 
Multiparticle Collision Description of the Solvent 
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A mesoscopic coarse-grain model for computationally-efficient simulations of biomembranes is 
presented. It combines molecular dynamics simulations for the lipids, modeled as elastic chains 
of beads, with multiparticle collision dynamics for the solvent. Self-assembly of a membrane from 
a uniform mixture of lipids is observed. Simulations at different temperatures demonstrate that 
it reproduces the gel and liquid phases of lipid bilayers. Investigations of lipid diffusion in differ- 
ent phases reveals a crossover from subdiffusion to normal diffusion at long times. Macroscopic 
membrane properties, such as stretching and bending elastic moduli, are determined directly from 
the mesoscopic simulations. Velocity correlation functions for membrane flows are determined and 
analyzed. 
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I. INTRODUCTION 



Biological membranes, formed by lipid bilayers, play 
a fundamental role in the function of biological cells 
and the theoretical description of their structure, proper- 
ties and dynamics is an important and challenging prob- 
lentii^. While powerful analytical theories exist for such 
systems2,ii, they treat membranes as mathematical sur- 
faces and do not reproduce the physical bilayer structure. 
Therefore, they are applicable only on scales that are 
much larger than the actual membrane thickness. More- 
over, they include phenomenological parameters which 
still need to be determined either from experiments or 
from microscopic simulations. For biological processes in- 
side a cell, the relevant length scales lie in the nanometer 
and submicrometer ranges. In order to consistently de- 
scribe processes involving micro- vesicles, membrane pro- 
teins and ion channels, theory and simulation that ac- 
count for the lipid structure of a bilayer are required. 

All-atom molecular dynamics (MD) simulations of 
lipid bilayers have been performed (see, e.g., Ref. [1- 
Q). However, they are costly and limited to relatively 
small systems and short time scales. While short-time- 
scale simulations are sufficient for the exploration of some 
aspects of membrane dynamics, there are many impor- 
tant biochemical processes which occur on longer time 
scales. For instance, it is known that characteristic 
mechanochemical motions in proteins, essential for their 
enzyme or motor functions, usually require milliseconds 
or more for their completions. Hence, microscopic inves- 
tigations of biomembranes with protein inclusions are be- 
yond the capacity of all-atom MD simulations. Further- 
more, such simulations are also too slow to microscop- 
ically reproduce the self-assembly of vesicles, structural 
instabilities of membranes or the effects of slow hydrody- 
namical modes on membrane dynamics. 



This has prompted the development of a variety 
of coarse-grain simulation methods for biomembranes, 
which are still able to resolve important aspects of the 
lipid bilayer structur&iSr— . Typically, a lipid molecule is 
modeled as a chain comprising one hydrophilic and sev- 
eral hydrophobic beads connected by elastic springs; each 
of these beads corresponds to a certain atomic group. In 
coarse-grain solvent descriptions, the solvent molecules 
are also represented by groups of atoms. 

In implicit solvent modelsi^r— , the solvent particles 
are not actually included in a simulation and hydropho- 
bic effects due to the presence of such particles are taken 
into account through the use of a tunable interaction po- 
tential between the lipids. Such a simplification results in 
a computational speed-up, making simulations of large- 
scale membrane instabilities possible^. However, in such 
solvent-free models the coupling of biomembranes to hy- 
drodynamic flows, as well as the hydrodynamic interac- 
tions mediated by the solvent, cannot be described. 

In explicit solvent models employing dissipative parti- 
cle dynamics (DPD), the solvent particles are included 
into the dynamical description, but actual molecular in- 
teractions between them are replaced by effective soft 
interaction potentials, so that the particles are allowed 
to penetrate one another. The use of a soft-core poten- 
tial for the solvent and lipids makes it possible to employ 
much larger molecular dynamics integration time steps 
compared with those in all-atom MD simulations; there- 
fore, substantially accelerating the computatio n^"^ i 
Nonetheless, further acceleration is desirable. 

A major portion of the computational time in ex- 
plicit solvent models is spent simulating the dynamics 
of the large number of solvent molecules in the system. 
This suggests that it is desirable to construct a coarse- 
grain dynamical scheme that treats the solvent part of 
the dynamics efficiently. Such a scheme is provided by 
multiparticle collision (MFC) dynamics^^i^i. In this ap- 
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proach, solvent particles, representing coarse-grained real 
molecules, free stream and undergo effective multiparti- 
cle collisions at discrete time moments. The collision 
and streaming rules are formulated in such a way that 
the mass and momentum conservation laws are satisfied. 
These rules can be constructed so that the dynamics is 
either micro-canonical and preserves the phase-space vol- 
ume or is canonical at constant temperature. MPC dy- 
namics has been applied to a variety of problems where 
fluid micro-flows were essential and there were interac- 
tions between fluids and macromolecules. Reviews of this 
method are availabl o^^i^^ . 

MPC dynamics has already been used for simulations 
of biomembranes. This method has been employed'^^ to 
study a micron-size vesicle under shear flow. In this work, 
the membrane was modeled as a triangulated surface de- 
scribed by vertices connected by tethers; the lipid bilayer 
structure was not resolved. In another study, a special 
color-collision rule was used to account for the interac- 
tion between the MCP solvent and the coarse-grained 
lipids'^^. In our investigation a coarse-grain description 
of the lipid bilayer, resolving membrane structure, is com- 
bined with MPC dynamics for the solvent. Interactions 
between lipids and solvent particles are explicitly taken 
into account. 

In Sec. ini the detailed formulation of the simulation 
method is given. Simulations for membranes at three dif- 
ferent temperatures are presented in Sec. IIIIl The simu- 
lations can reproduce a gel phase at low temperature and 
a liquid phase at higher temperatures. Density profiles 
for lipid particles across the membrane, lipid chain or- 
der parameters, and radial distribution functions of lipid 
head particles are determined and discussed. Through 
direct simulations, intra-membrane diffusion is explored 
and a subdiffusion regime on relatively short time scales 
is observed. In the next sections, our investigations fo- 
cus on the membrane in the liquid phase, important for 
biological applications. In Sec. IIVI self-assembly of a 
membrane from an initially uniform mixture of lipids is 
demonstrated. The surface tension coefficient is deter- 
mined from simulations on membranes of different sizes 
in Sec.|Vl By constructing and analyzing the power spec- 
trum of membrane height fluctuations, the elastic bend- 
ing modulus of the membrane is found, fluctuations of 
the membrane flow velocity are considered and velocity- 
velocity correlation functions are analyzed. The paper 
ends with conclusions and a discussion of the results. 



II. MESOSCOPIC MODEL FOR LIPID 
BILAYER DYNAMICS 

In this section we describe the mesoscopic coarse-grain 
model for the structure and dynamics of a lipid bilayer 
membrane in a solvent. The mesoscopic model uses a 
coarse-grain description of a lipid molecule as a collec- 
tion of linked molecular groups termed beads. In addi- 
tion, the solvent in which the lipids reside is treated at 




FIG. 1. The lipid chain (left) and its schematic representation 
as a rod (right). The lipid consists of four beads linked by 
elastic FENE bonds (solid lines) and straightened by elastic 
bonds (dashed lines). The first bead (dark gray/blue) is hy- 
drophilic. Three other beads are hydrophobic, the terminal 
hydrophobic bead is shown as light gray /blue 



a particle-based level where each effective point solvent 
particle represents a collection of real solvent molecules. 
The coarse-grained lipid molecules interact through inter- 
molecular potentials. The solvent particles also interact 
with the lipid beads through intermolecular potentials; 
however, the solvent particles interact among themselves 
through multiparticle collisions. There are no intermolec- 
ular interactions among solvent particles. The dynami- 
cal evolution of the entire systems, lipids plus solvent, is 
described by a hybrid dynamical scheme that combines 
molecular dynamics for all interacting particles with mul- 
tiparticle collision dynamics for the solvent. The fact that 
there are no explicit solvent-solvent molecule interactions 
is responsible for the computational efficiency of this dy- 
namical scheme. Below we provide a detailed description 
of the mesoscopic MD-MPC dynamical bilayer model. 

A. Lipid interactions 

A lipid chain comprises a hydrophilic head and a hy- 
drophobic tail. In common with many other coarse-grain 
descriptions, a lipid molecule is modeled as a set of beads. 
In our investigation, we adopt a four-bead representation 
of the lipid where the hydrophobic head (h) is modeled as 
a single bead and the hydrophobic tail (t) as three beads 
(see Fig. [Ija)). Below, we specify the interactions be- 
tween the beads in a lipid and between the lipids. These 
lipid interaction potentials have the same forms as in 
Cooke, et ali^. 

The interaction between two lipid beads is described 
by the truncated Lennard- Jones (LJ) potential. 
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where 6{r) is the Heaviside function and rij = \ri — rj\ 
is the distance between the beads i and j. The cutoff 
length = 2i/6cr is chosen in such a way that there is a 
short-distance repulsion but the long-distance attraction 
is absent. The strength of the interaction between beads 
i and j takes the value eaa', where a, a' £ {h,t} if bead 
i is of type a and bead j is of type a' . 

Two neighboring particles in a lipid chain are linked 
by a FENE bond'^-, described by the potential 
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where ri is the distance between the beads and roo = 
1.5 a is the maximum distance allowed by the FENE 
bond. In all simulations, we have chosen the spring con- 
stant as kbond = 20e/tT^, so that at equilibrium the length 
of the FENE bond is close to a. Bending rigidity of a 
lipid chain is modeled by introducing additional springs 
connecting next-nearest neighbor beads and is described 
by the bending potential, 
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Vbend = -zhend{r2 " 4cr) 



(3) 



where is the distance between such two beads. The 
spring constant is kbend = 2.5e/ cr^ and the natural length 
is 4 a. For a slightly bent lipid chain with FENE bond 
length (7, this potential reduces to ^kbendC^S^ , where tt — 
9 is the angle between two neighboring FENE bonds. 
Hence, it provides a bending stiffness of 2.5e to the lipid. 

Hydrophobic effects are responsible for the aggregation 
of lipids into a membrane. These were taken into account 
by adding an attractive potential between beads that be- 
long to different lipid tails. The effective interaction was 
chosen to be 
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where Waa' = Wtt and taa' = ^tt when interactions be- 
tween two tail beads from different lipid chains are con- 
sidered. This lipid model was originally constructed to 
describe a lipid membrane in the absence of solvent^. 
Since our simulation contains explicit, albeit effective 
point solvent molecules, the parameters that enter this 
model were altered (see below) to account for the explicit 
presence of the solvent molecules. 



B. Lipid-solvent interactions 

The solvent particles interact with the lipid beads 
through intermolecular potentials. The interaction be- 
tween a solvent particle and a lipid tail bead is also given 



by Eq. ([T]) with the same cutoff length Tc, but with a 
different interaction strength Caa' = Est- This interac- 
tion is purely repulsive; it accounts for hydrophobic ef- 
fects. The interaction between a solvent particle and a 
lipid head bead is given by Eq. (j4]) with Waa' = Wsh and 
^aa' = Cs/i- This interaction is repulsive at r < Tc and at- 
tractive for Tc < r < rc + w^h , so that hydrophilic effects 
are taken into account. 



C. MD-MPC dynamics 

The system consists of Nl lipid molecules and Ns sol- 
vent molecules. Since there are no explicit solvent-solvent 
interactions, the total potential energy of the system, Vt, 
may be written as the sum of interactions within the 
single lipid molecules, Ve, interactions among different 
lipid molecules, Vu, and lipid-solvent interactions. Vis- 
Vt — Vi+Vee+Vis- Instead of explicit interactions among 
solvent molecules, their interactions are treated by multi- 
particle collision dynamics^^. Hybrid MD-MPC dynam- 
ics combines molecular dynamics segments of evolution 
with effective multiparticle solvent collisions at discrete 
time intervals r to obtain the time evolution of the entire 
system in the following way2i: 

Given that the total potential energy of the entire sys- 
tem is Vt, Newton's equations of motion are used to 
evolve all particles for a time interval t. Note that be- 
cause there are no solvent-solvent interactions this MD 
trajectory segment can be simulated efficiently, even for 
large systems containing many solvent particles. At time 
T multiparticle collisions among solvent molecules take 
place. To carry out such collisions, the solvent particles 
are sorted into the cells of a simple cubic lattice and par- 
ticles in the same cell exchange momentum with each 
other while the total momentum in the cell is conserved. 
We employ the constant temperature version of MFC dy- 
namics^-. If the mean velocity of the solvent particles in 
the cell ^ is V^, the collision event of the i-th particle 
inside this cell is modeled by updating its velocity, V;, so 
that the new velocity, v^, is given by 
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where the components of v[°" are chosen as Gaussian 
random numbers with zero mean and variance fcsT/m, 
is the number of solvent particles in the cell ^ and 
the summation is performed over all solvent particles in 
this cell. Since the mean free path of the solvent parti- 
cles in our simulation was small compared with the size 
of a MFC cell, we used random grid-shifting'^^'~ to im- 
plement the MFC step. This sequence of MD and MFC 
steps is repeated to evolve the entire system. The proper- 
ties of such MFC dynamics have been discussed in detail 
in reviews where further applications can be found™—. 



4 



D. Simulation details 

The characteristic interaction energies between differ- 
ent types of beads were eht = 1 e, Qih — ^tt — 0.5 e, 
Esh = 0.05 e and est = 2.0 e. The attraction ranges for 
tail-tail and solvent-head interactions were chosen such 
that Tc -\- Wtt = 2.6 a and Tc -I- Wsh = 1-65 a. All particles 
and beads had equal mass m. The simulations were car- 
ried out in a cubic box of size 25 a x 25 a x 25 cr with peri- 
odic boundary conditions. The lateral size of a MPC cell 
was uq = (T. The system contained 1000 lipid chains and 
56624 solvent particles. On average, the solvent number 
density in the bulk was equal to five. 

The initial velocities of all particles were Gaussian dis- 
tributed with zero mean and variance fcsT/rn for each 
component. For the MD trajectory segments, Newton's 
equations of motion were integrated using the velocity- 
Verlet algorithn t^ wit h a time step of St — 0.005 to, 
where to = -y/mcr^/e, and the MPC time step was 
r = 0.2 = 40 6t. The initial configuration of the mem- 
brane was prepared by arranging the lipids as a bilayer 
in the xy-plane, with the hydrophilic particles facing the 
outer surfaces while the solvent particles were randomly 
distributed in the rest of the simulation box. Simulation 
data was gathered after the system had evolved for 10^ St, 
so that thermal equilibrium was established. Depending 
on the physical quantity under investigation, time aver- 
ages were taken over time intervals up to 10^ St. 

Results will be reported below in dimensionless simula- 
tion units except where connections with physical length 
and time scales are made. We have chosen a to be the 
unit of length and m the unit of mass. The characteristic 
interaction energy between a lipid head and a lipid tail 
bead, eht = e, was taken to be the unit of energy. Time 
will be reported in units of St. 



III. MEMBRANE PROPERTIES AT 
DIFFERENT TEMPERATURES 

Self-assembled lipid membranes are known to have a 
rich phase behavior^^-. At higher temperatures, the lipids 
in the membrane are not ordered and the membrane is in 
the so-called liquid phase. As the temperature decreases, 
the membrane undergoes a transition to a gel phase in 
which the lipid chains show nematic order. In simula- 
tions at three different temperatures, we observed various 
bilayer structures, which were analyzed by determining 
the vertical density profile, the lateral radial distribution 
function, the chain order parameters, and the in-plane 
diffusion constant of the lipids. 

Examples of membrane structures observed in our sim- 
ulations are shown in Fig. [51 We have chosen to visualize 
the lipids using the rod representation shown in Fig.[ljb). 
The FENE bonds are displayed as gray solid rods, with 
only the hydrophilic head beads (dark blue) and the ter- 
minal bead of the hydrophobic tail (light blue) explicitly 
shown. 




FIG. 2. Membrane structures at three different temperatures: 
(a) ksT/e = 0.4, (b) ksT/e = 1.0 and (c) kaT/e = 2.0. 
The rod representation is used to display lipids. The solvent 
particles are not shown. 

At ksT/e = 0.4 (Fig. [IJa)), the lipids are mostly 
straight and relatively well ordered. Two domains can 
be seen in the figure. In the majority domain, the lipids 
are tilted and roughly parallel to one another. In the 
smaller domain in the right part of the membrane, the 
lipids are less ordered and the structure is similar to that 
seen at the higher temperatures. At ksT/e = 1.0, the 
tilted ordered structure is not observed and the orien- 
tation of lipids is less ordered (Fig. [DJb)). Nonetheless, 
there is still a well-defined midplane which separates two 
lipid monolayers. Moreover, it is clear that, on average, 
the lipid chains are perpendicular to the midplane of the 
membrane. When the temperature is further increased 
to ksT/e — 2.0, an irregular structure is found where 
the orientational order of the lipid chains is weak and the 
lipid head particles penetrate into the membrane interior, 
so that the bilayer midplane and the interface between 
the lipids and the solvent are less defined (Fig. HKc)). 

To quantitatively characterize the bilayer structures at 
various temperatures, the orientational order parameter, 
the in-plane radial distribution function and the vertical 
density distributions of lipid chains can be used. 

The orientational order parameter is defined as = 
^(3 cos^ 6g — l) where the bracket (...) denotes a canonical 
equilibrium average. For the i-th lipid chain, cos de = 
Ti-n, where ff is the unit vector pointing from the last tail 
bead to the lipid head and fi is either the unit normal to 



the upper or lower monolayers. The orientational order 
of the chains decreases as S diminishes. When = 1, 
all lipid chains are aligned parallel to the bilayer normal. 
On the other hand, S — implies that, on average, there 
is no correlation between the directions of the lipids and 
the bilayer normal. 

In our simulations, the orientational order parameter 
was determined by averaging over all lipids and over 
1000 bilayer configurations separated by 200 St. We 
found that S = 0.54 for the membrane at tempera- 
ture kBT/e = 1.0, typical for a membrane in the liquid 
phaseSS. At the higher temperature ksT/e = 2.0, the 
order parameter drops to = 0.21, thus indicating a 
more disordered orientational structure. One might have 
expected that the lipid chains would have been more or- 
dered at the lower temperature fc^T/e = 0.4. However, 
the orientational order parameter actually decreases to 
S = 0.46, since most of the chains are then tilted and 
therefore their direction deviates from the bilayer nor- 
mal. 

To better characterize chain orientational order at the 
temperature fc^T/e = 0.4, we have chosen a domain 
where the lipids were tilted and introduced the unit vec- 
tor rit pointing along the average direction of the tilted 
lipids. In this domain, cos 9i was determined by comput- 
ing the inner product of fit and the unit vector of the 
lipid. When cos 9£ was defined in this way, we found that 
S = 0.96, confirming that the orientational order of the 
membrane was even higher at this lower temperature. 

We have also determined the in-plane radial distribu- 
tion function .g(r||) of lipid head beads, 



where /9(r|||0) is the average two-dimensional density 
of head beads at a projected distance r|| on the xy- 
plane from a given head bead and p is the average two- 
dimensional density of lipid-head beads. To compute this 
property averages were taken over all lipid head beads in 
500 bilayer configurations separated by 2000 6t. Figure [3] 
shows radial distribution functions g{r\\) at three differ- 
ent temperatures. When ksT/e = 0.4, the radial distri- 
bution function has several peaks extending to r|| w 5 
and the separation between the peaks is close to the size 
of a lipid bead. At fcsT/e — 2.0, the fact that the ra- 
dial distribution function is not vanishing when r|| < 1 
suggests that the positions of two lipid head beads, pro- 
jected on the xy-plane, overlap due to the presence of 
lipid head beads in the interior of the membrane. As 
temperature increases, the in-plane correlations become 
weaker indicating that the membrane is less structured 
at higher temperatures. 

Another important statistical property is the vertical 
density profile of lipid particles. Figure S] displays a cut 
through the simulation box showing the vertical struc- 
ture of the membrane and surrounding solvent particles. 
To determine the vertical profiles, the simulation box was 
divided into 250 slices in the z-direction; each slice had 




FIG. 3. Radial distribution functions of lipid head beads in 
the membrane at three different temperatures: (a) ksT/e — 
0.4 (dashed black line), (b) ksT/e = 1.0 (thick blue line) and 
(c) ksT/e = 2.0 (thin red line). 




FIG. 4. A cut through the simulation box showing the vertical 
structure of the bilayer and solvent particles at ksT/e = 1.0. 



thickness 0.1. The time-averaged density profiles of sol- 
vent {ps)^ lipid head {ph) and lipid tail {pt) beads for each 
slice at different temperatures were computed, with the 
average over all system configurations up to 10^ 5t after 
the system reached equilibrium. The results are shown 
in Fig. [5l 

When ksT/e = 0.4, the density profile of lipid tail 
beads consists of several sharp peaks, each of which corre- 
sponds to the vertical position of one lipid-tail bead, thus 
indicating a well-ordered vertical arrangement for the 
beads along a chain and small membrane shape fluctua- 
tions (Fig. [SJa)). As temperature increases, a smoother 
profile for the lipid tail density is observed, showing that 
the beads along a lipid chain are less ordered and thermal 
fluctuations of the membrane shape are more significant 
(Fig-EIb)). At ksT/e = 2.0, one can see that the density 
of lipid head beads in the interior of the bilayer becomes 
significant. Moreover, the distribution of lipid tails is 
also broader and a larger overlap with the distribution 
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FIG. 5. Vertical density profiles for hydrophilic head beads 
(ph, thick black line), hydrophobic tail beads {pt, dashed blue 
line) and solvent particles {ps , thin red line) at three different 
temperatures : (a) ksT/e = 0.4, (b) ksT/e = 1.0 and (c) 
ksT/e = 2.0. The scale is different for the solvent density 
profile. 

of lipid- head beads is observed (Fig. EKc)). These data 
again indicate a more disordered bilayer structure, close 
to the onset of membrane dissociation. 

The vertical density profile of the lipid beads, the ra- 
dial distribution function of the lipid-head beads and the 
lipid chain order parameter give us information on the 
equilibrium organization of the membrane. As the equi- 
librium structure of the membrane changes with temper- 
ature, the dynamics of individual lipids should also be 
affected. 

To investigate lipid diffusion, we computed the in- 
plane mean square displacement of the lipids, MSD(t) = 
([r^ii (t)— r£ II (0)]^), where || is the position of the center 
of mass of a lipid projected on the xy-plane. The aver- 
ages were taken over all lipids at every 1000 St and the 
MSD was computed from system trajectories of length 
up to 10^ 6t. Depending on the time domain, both dif- 
fusive and subdiffusive types of behavior of lipids were 
found. 

For all three temperatures lipid diffusive motion was al- 
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FIG. 6. Diffusion of lipids in the membrane. Log-log plots of 
the mean square displacements (MSD) of the center of mass 
of a lipid are shown as functions of time for (a) fcsT/e = 0.4, 
(b) ksT/e = 1.0 and (c) fcsT/e = 2.0. The dashed and solid 
straight lines are linear fits for the subdiffusive and normal 
diffusive regimes, respectively. 



ways observed in the long-time regime, so that MSD(t) — 
4Dt, where D is the diffusion constant. At ksT /e — 0.4 
(Fig. IHa)), the diffusion constant is D = 7.87 x 10"^, 
implying that a lipid moved over a distance approxi- 
mately equal to the size of a coarse-grained lipid-head 
bead within the simulation time of t ~ IQ^St. At a higher 
temperature ksT/e ~ 1.0 (Fig. [BKb)), the diffusion con- 
stant is 13 = 7.86 X 10~^, which is about 10 times larger 
than that at the lower temperature. The diffusion con- 
stant further increases to D = 1.86 x 10"'' at /csT/e = 2.0 
(Fig. Efc)). Thus, lipid diffusion in the bilayer depends 
strongly on temperature. 

In addition to normal lipid diffusive dynamics in the 
long-time limit, subdiffusive motion was found at inter- 
mediate times, so that MSD(t) ^ t", where a < 1 is the 
subdiffusive exponent. At /csT/e = 0.4 (Fig. [HUa)), we 
found a = 0.44 over the time up to the crossover time 
tc ~ 3900 St. As the temperature increases, the subdiffu- 
sion exponent a grows to 0.82 and 0.9 for ksT/e =1.0 
and ksT/e = 2.0, respectively. Similar subdiffusive be- 
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havior has been recently observed in all-atom MD simu- 
lations of lipid membranesSii, where in-plane motions of 
lipids were subdifFusive on the time scale of nanoseconds. 

Our simulation data suggests that the membrane at 
/csT/e = 1.0 is in liquid phase, with a well-defined 
solvent-lipid interface and bilayer midplane. Since liquid- 
like membrane states are most relevant for real biological 
membranes, all our subsequent simulations reported be- 
low were carried out at k^T je = 1.0. 

At this temperature, the mean vertical distance be- 
tween the peaks in the lipid head density profiles for two 
monolayers (see Fig. EJb)), which can be chosen as the 
mean thickness of the membrane, is close to 5.5 <j. Com- 
paring this to the thickness of a typical real membrane 
(between 4 and 6 nm) , we can identify the unit length in 
our simulations to be cr ~ 1 nm. The time unit Et can be 
roughly estimated by comparing the typical lipid lateral 
diffusion constant in experiment^ (« 4 /im^/s) with our 
computed values. We find 8t ~ 20 ps. Hence, the phys- 
ical size of the simulation box is about 25 nm and the 
total simulation time is about 20 /iS for the simulation 
with 10^ <5t. 



IV. SELF-ASSEMBLY OF THE MEMBRANE 

To verify that the lipid bilayer is indeed thermodynam- 
ically stable at /csT/e = 1.0, simulations of the mem- 
brane self-assembly process at that temperature were 
performed. The initial condition was taken to be a ran- 
dom mixture of lipids and solvent particles. To prepare 
this initial configuration, the following procedure was em- 
ployed: starting with an equilibrium lipid bilayer, the 
attractive potential interactions between the lipid tail 
beads were switched off. In addition, attractive interac- 
tion potentials between the tail beads and solvent parti- 
cles, with est = 0.1 and -I- Wst = 2.6, were introduced, 
making lipids more hydrophilic. After 2 x 10^ (5t, the 
lipid and solvent particles were found to be uniformly 
distributed within the simulation box. 

Starting from this uniform initial configuration, sim- 
ulations with the full potential model were performed. 
FigurejTIJa) shows the initial configuration at f = where 
the lipid chains are uniformly distributed inside the sim- 
ulation box. As time evolves, the lipids quickly aggre- 
gate forming small segments (Fig. El^b)). These small 
segments gradually merge into large branched clusters 
(Fig. El^c)). Then, a slow rearrangement process takes 
place leading to a single bilayer structure with a large 
hole in the center (FiglTJd)). Later, this large hole slowly 
shrinks to a small pore (Fig[7l^e)) and eventually closes 
after 3 x 10^ bt (FigE^f)). The formation of the lipid 
bilayer and the close-up of the membrane pore suggest 
that the uniform, flat membrane is indeed a thermody- 
namically stable structure at fc^T/e = 1.0. (See on-line 
video of the evolution process.) 




FIG. 7. Self-assembly of the lipid bilayer. Six configurations 
at time moments (a) (5t, (b) 6000 &t, (c) 18000 5t, (d) 160000 
(5t, (e) 246400 &t and (f) 300000 &t are shown. The initial 
state (a) corresponds to the uniform mixture of hpids and the 
solvent. (Enhanced online) 



V. MACROSCOPIC MEMBRANE PROPERTIES 
AND VELOCITY CORRELATION FUNCTIONS 

Macroscopic theories of biomembranes are formulated 
in terms of elastic deformable surfaces. There are 
also hydrodynamic descriptions which treat the mem- 
branes as two-dimensional fluids immersed in the three- 
dimensional solvent. Below we show that our simulation 
data is consistent with such macroscopic theories. More- 
over, our mesoscopic simulations allow us to determine 
some of the characteristic properties of the lipid bilayer 
membrane. The analysis is restricted to the liquid-phase 
membrane at ksT/e — 1.0. 

The surface tension of the membrane can be obtained 
by considering the membrane stretching energy. In 
the regime where Hookian elasticity theory holds, the 
stretching elastic energy of a membrane with area A is 
given by 



where Ka is the stretching modulus and Aq is the refer- 
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FIG. 8. Dependence of the surface tension 7 on the membrane 
area A. The first five data points were used to determine the 
membrane stretching modulus. 



and the area Aq corresponding to the tensionless mem- 
brane. We found that Aq ~ 567 cr^ and — 33.9 e/cr^, 
In Sec. mil we noted that c ~ 1 nm. Since our simula- 
tions were performed at ksT / e = 1.0, we have e — ksT. 
Therefore, the computed stretching modulus is approxi- 
mately Ka = 33.9 fc^r/nm^. This is comparable to the 
values observed for typical liquid-like membranes^, i.e. 
if^ = 50 - 70 ksT/mr?. 

In the macroscopic continuous approact^, the Helfrich 
free energy of the membrane is 



^ JLxL 



(11) 



ence area corresponding to a tensionless membrane. The 
membrane surface tension is 7 = dEg/dA and therefore 
the equation 



1 = Ka 



A-Aq 
Aq 



(8) 



holds. In our simulations, we determined the equilib- 
rium surface tension of the membranes with different ar- 
eas confined to boxes of different lateral sizes. From these 
measurements, the stretching modulus could also be ob- 
tained. 

The surface tension of a membrane was determined 
in our simulations from the pressure tensor using the 
relation^, 



P. 



1 



(9) 



where the bracket denotes an equilibrium canonical aver- 
age and Lz is the linear size of the simulation box in the 
z-direction. The diagonal elements of the pressure tensor 
are defined as 



P -1 
~ V 



(10) 

where V is the volume of the simulation box. The sum- 
mations in this equation are taken over all particles, in- 
cluding both the solvent and the lipids. The z^-component 
of the distance between two particles, i and j, is rij^i, and 
the force acting between them is When both par- 

ticles are solvent, fij^^, = 0; otherwise, fij^i, is evaluated 
through the actual interaction potentials. The equilib- 
rium average in Eq. Q was computed by a time average 
over an interval of 10^ 5t. 

The computed values of the surface tension for dif- 
ferent membrane areas are shown in Fig. |8l When the 
surface tension is small, it depends approximately lin- 
early on the area A. By fitting this linear dependence 
to Eq. ([5]), we determined the stretching modulus Ka 



where /i(r) = hqC^"^'^ is the local height of the mem- 
brane measured with respect to the reference plane. As 
implied by the energy equipartition theorem, the power 
spectrum of membrane height fluctuations should there- 
fore be 



S{q)^L^\hl\) 



Kq + jq 



(12) 



where n is the membrane bending modulus and 7 is 
again the surface tension. There exists a characteris- 
tic wavenumber qc = \/jJk separating two different 
regimes. When q <^ q^ the power spectrum is S{q) ^ q^^ 
and the dominant contribution comes from the mem- 
brane tension. For q ^ Qc, the power spectrum is 
S{q) ~ and the dominant role is played by the bend- 
ing elasticity. 

To determine membrane height fluctuations, a bi- 
layer configuration from a simulation was taken at every 
5000 so that, in total, 1000 bilayer configurations were 
recorded. In each bilayer configuration, the membrane 
was divided into a grid of 10 x 10 cells. The membrane 
height of each cell was further determined by taking the 
average of the positions of end beads in the hydrophobic 
lipid tails. In this way, local heights could be determined 
at all grid points. Performing a fast Fourier transform for 
the membrane heights, the power spectrum S{q) could 
be determined for each bilayer configuration. By averag- 
ing over all 1000 bilayer configurations, the mean power 
spectrum was obtained. Note that, based on our simu- 
lations, the power spectrum could only be computed in 
the range qmax < q < qmin- Here qmax = Stt/Z is twice 
the linear size of a grid cell, close to the bilayer thickness, 
and qmin = Stt/X, where L is the linear dimension of the 
simulation box. In our simulations, we had L — 25(7 and 
I = 5 a, so that qmax — 1-25 and qmin — 0.25 a^^. 

Figure [9] displays the numerically determined power 
spectrum. The solid line shown in Fig. [5] is obtained by 
least-squares fitting using Eq. (|12p with the membrane 
tension value 7 = 3.21 fc^T/nm^ taken from the constant 
surface area simulations (Fig.[8|). As a result of data fit- 
ting, the membrane bending modulus was found to be 
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1 

FIG. 9. Power spectrum S{q) of membrane height fluctua- 
tions. The solid line is the best fit of the simulation data, 
using the theoretical dependence (Eq. (fT^ ). 

K ~ 12.56 e ~ 12.56 ksT. Using this value of the bend- 
ing modulus K, and the previously determined value of 
the surface tension 7 for the membrane, the character- 
istic wavenumber = 0.5 a^^ could be obtained. This 
wave number lies in the middle of the computed power 
spectrum, indicating that our simulations are able to re- 
produce both the tension-dominated and the bending- 
dominated regimes. Typical experimental values of the 
bending modulus k for lipid membranes lie^- between 10 
and 20 ksT. Hence, we can again notice that the mem- 
branes in our simulations are similar in their physical 
properties to real biological membranes. 

Finally, we consider flow dynamics of lipids in the 
membrane. In the classical study by Saffman and 
Delbriick'*'^, the membrane was treated as a two- 
dimensional (2D) simple fluid embedded in a three- 
dimensional (3D) solvent. When a lipid moves in the 
membrane, its momentum may be transferred not only 
to the neighboring lipids, but also to the solvent. How- 
ever, estimates show^^iii that, on length scales shorter 
than a micrometer, hydrodynamic coupling between the 
membrane and the solvent is not significant and, on such 
scales, the membrane can be approximately treated as a 
2D fluid. 

The longitudinal and transverse velocity correlation 
functions of lipid flows are 

CL{x,t) = {va:{xo,ya,to)va;{xQ + x,yo,to + t))^^^_^^^^^, 
CT{x,t) = {vy{xQ,y(),ta)vy{xo + x,yo,to+t))^^^ 

' as) 

where angular bracket {■ ■ ■)xo,yo,to denotes an average 
over the positions xq, yo, time to and realizations. The 
hydrodynamic velocity field •v{x,y,t) is defined by tak- 
ing the average of the instantaneous velocities of all 
lipids within a certain membrane area element. In our 
simulations, the membrane was divided into a grid of 
10 X 10 of cells and the hydrodynamic velocities were ob- 
tained by averaging the in-plane lipid velocities in each 
cell. The products v.j;{xQ,yQ,tQ)vx{xo + x,yo,to + t) and 
Vy(xQ, yo, tQ)vy{xo + X, yo, to + t) were determined for all 




FIG. 10. Time dependence of the longitudinal (a) and trans- 
verse (b) velocity correlation functions for three different sep- 
arations: a: = (solid lines), x = 0.25 (dashed lines), and 
a; = 5 (dotted lines). 

grid points (a;o,yo) at every MD step, and the correla- 
tion functions CL{x,t) and CT{x,t) were computed by 
taking the average of these products over all grid points 
(xo,2/o) and over 1000 St. Subsequently, the results were 
additionally averaged over an ensemble of 20 independent 
realizations. 

Figure [10] shows the dependences of CL{x,t) and 
Cxixji) on time for three different values of the dis- 
tance X. The peak in Cl (2.5, t) is found 40 St later than 
the peak in (7^(0,^), suggesting that it takes 40 St for a 
fluctuation of velocity Vx to be transported over a dis- 
tance X — 2.5. Similarly, it takes about 25 St for the 
fluctuations of Vy to be transported in the x-direction 
over such distances. This is much faster than the time, 
x^/D ~ 8 X 10'^ St for x = 2.5, needed for the lipids to 
diffuse over the same distance. Therefore, we conclude 
that velocity fluctuations are transported by collective 
lipid flows, not by the diffusion of single lipids. 

We can also consider the time integrals of the velocity 
correlation functions, 

/■OO 

GlIx) = / CL{x,t) dt, 

(14) 

Gt{x) = / CT{x,t) dt. 
Jo 

They are determined by the pair mobility tensor which 
describes the velocity response of one fluid element due 
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FIG. 11. Longitudinal (circles) and transverse (squares) cor- 
relation functions Gl{x) and Griy)- The solid and dashed 
lines show the respective logarithm approximations given by 
Eq. 



to the motion of another element in the fluid^ . Such 
responses are given by the Green function of the Stokes 
equation. The behavior of the Green functions depends 
on the dimensionality of the fluid. For three-dimensional 
fluids, the functions fall as 1/r with the distance r. In 
contrast to this, logarithmic distance dependence is char- 
acteristic for two-dimensional fluids. 

As suggested by Saffman and Delbriick^, biomem- 
branes can be viewed as 2D fluids of lipids which are 
immersed in a 3D solvent. On length scales typical for 
our simulations, viscous coupling between the membrane 
and the solvent is negligible. Assuming that the mem- 
brane is a planar 2D fluid, expressions for the longitudinal 
and transverse velocity fluctuations can be derived from 
the pair mobility tensor'*'*. Thus, one gets 

Gl{x) = -CHx/R,), 
GT{x)^-C[l + \n{x/R,)], 

where C is a constant prefactor and Rc is a cutoff length 
which is typically on the micron scale. These approxi- 
mate expressions hold for distances x < Rc- On longer 
length scales, momentum diffusion into the bulk solvent 
becomes significant and a crossover to the behavior char- 
acteristic for 3D systems should take place. Note that 
for a finite system, Rc should be approximately equal to 
the linear system size^. 

Figure [TT] displays Gl{x) and Gt{x), the longitudinal 
and transverse correlation functions, determined in our 
simulations. The solid and dashed lines show best fits 
using the logarithmic approximations (jlSp with C ~ 1.6 x 
10"'^ and Rc — 20. Good agreement is found indicating 
that the lipid flows in our simulations were indeed well 
described in terms of 2D hydrodynamics and that the 
leakage of lipid momentum into the solvent was negligible 
on the length scale of our system. 



VI. DISCUSSION AND CONCLUSIONS 

We have presented and tested a coarse-grain simulation 
method for biomembranes. In common with other coarse- 
grain methods, individual lipids were modeled as short 
chains of beads linked by elastic bonds, and the solvent 
was explicitly included using multiparticle collision dy- 
namics. Our method differs from other investigations2fii^ 
of lipid membrane where MFC dynamics for the solvent 
was employed in that we account both for the structure of 
lipid bilayer and include explicit lipid-solvent hydropho- 
bic and hydrophilic interactions. 

The interaction parameters of the model were chosen 
to reproduce the behavior of typical real lipid bilayers. 
Thus, we could follow in our simulations the self-assembly 
of a membrane starting from a uniform mixture of lipids 
and solvent. We could also reproduce various structural 
states of lipid bilayers at different temperatures, includ- 
ing the gel phase at the lower temperature and the liquid 
phase at the higher temperature. 

Statistical properties of collective modes of the liquid 
state of the membrane were studied. By varying the 
membrane area, the membrane surface tension was deter- 
mined and the lateral stretching modulus were obtained. 
The bending modulus of the membrane was then derived 
from the power spectrum of membrane height fluctua- 
tions. The results show that the elastic properties of our 
model membranes are comparable to those of a typical 
real lipid bilayer. 

Hydrodynamics of membrane flows was numerically 
investigated by computing correlation functions of the 
lipid velocity field. We found that the velocity fluctu- 
ations are not due to the diffusion of single lipids but 
are propagated by collective hydrodynamic modes. The 
computed velocity correlation functions show logarithmic 
spatial dependence, suggesting that, on the length scale 
of our simulations, the lipid bilayer could be considered 
as a 2D viscous fluid with little momentum diffusion into 
the bulk solvent. 

Our simulation method has a number of advantages. 
By modeling the solvent using multiparticle collision dy- 
namics, one does not need to expend computational 
power to calculate forces acting between solvent particles, 
as in MD and DFD simulations. Thus, the simulations 
could be substantially accelerated. 

Another important feature in our simulations was that 
the lipid-lipid and lipid-solvent interactions both con- 
tained short-range hardcore repulsion. Therefore, crowd- 
ing effects in the lipid membrane could be well repro- 
duced, as seen in the observed short-time subdiffusive 
motion of single lipid chains. This effect was previously 
reported in an all-atom MD study22,, but the long-time 
normal diffusive regime of single lipid chains was not 
found. 

Finally, we would like to point out that it is possible to 
combine our fast coarse-grain descriptions of membranes 
and solvent with coarse-grain simulations for proteins^. 
Such structurally-resolved numerical investigations of in- 
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dividual protein machines in biomembranes, as well as 
the collective dynamics of such protein machines, will be 
presented in future work. 
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